Statistics in Medicine

Selection bias is a common concern in causal inference. Two mechanisms can give rise to selection bias: collider bias or stratifying on an effect modifier of the exposure-outcome causal effect. Research into the impact of selection bias on instrumental variable (IV) analyses has focused primarily on collider bias. Here, we show that IV estimates can be biased if selection is an effect modifier (or caused by an effect modifier) of the exposure-outcome effect. In a simulation study, we investigate the size and direction of the bias when estimating effects on the mean difference or odds ratio scale. Bias can also occur due to heterogeneity in the instrument-exposure relationship, although its magnitude is typically smaller. Finally, we discuss methods to deal with this bias in practice, including inverse probability weighting (IPW).

The rapid growth of technologies such as wearable devices and high-throughput sequencing has generated unprecedented volumes of longitudinal data, often with multiple correlated features. This trend highlights the need for robust clustering methods tailored to multi-feature longitudinal datasets. Although previous reviews have benchmarked clustering methods in this context, most evaluations relied heavily on simulated and overly simplified datasets, which substantially limit their interpretability and applicability to real-world data. Missing data is also ubiquitous in longitudinal data, yet understudied in the context of clustering. To address these gaps, we systematically evaluate clustering methods across three distinct types of real-world longitudinal datasets: (1) a high-dimensional, low-timepoint metagenomic study, (2) a low-dimensional, high-timepoint multi-site repeated-measure study, and (3) a high-dimensional, high-timepoint repeated-measure study. We further investigate how missing data and different imputation strategies influence clustering performance. In addition, we introduce a novel resampling framework based on the empirical cumulative distribution function (eCDF) and a copula. This framework enables us to simulate multi-feature longitudinal data that preserves the correlation structure of real datasets, while also allowing for controlled variations such as sample size. Our simulation framework naturally accommodates missing data, thereby providing a flexible tool for benchmarking and comparing clustering methods in more realistic longitudinal settings.

AO_SCPLOWBSTRACTC_SCPLOWO_ST_ABSBackgroundC_ST_ABSLung donor-recipient matching requires balancing recipient medical urgency with expected post-transplant benefit. Conventional prediction models estimate outcomes for historical donor-recipient pairs but do not directly compare alternative candidates for the same donor. We developed an offline contextual-bandit framework to assess the feasibility of learning donor-centered lung-matching policies from historical registry data. MethodsWe used linked United Network for Organ Sharing (UNOS)/Organ Procurement and Transplantation Network (OPTN) Standard Transplant Analysis and Research files to identify 8,529 lung transplant allocation events from 2023 through 2025. Each event included the historical recipient and nine date-matched, ABO-compatible pseudo-candidates. The reward combined 72-hour respiratory-support-free survival (weight, 0.75) with a normalized initial Composite Allocation Score waitlist medical-urgency score (weight, 0.25). We trained a pessimistic neural lower-confidence-bound contextual-bandit policy (NeuraLCB) and evaluated it on held-out donor events using matched-action evaluation and self-normalized inverse propensity scoring (SNIPS). ResultsThe reconstructed dataset contained 85,290 donor-candidate rows. In the held-out set of 1,706 allocation events, the observed clinician reward was 0.534 (95% CI, 0.518-0.549). NeuraLCB agreed with the historical recipient in 174 events (10.2%), with a matched-action reward of 0.616. The SNIPS-estimated value of the NeuraLCB policy was 0.649 (95% CI, 0.568-0.717), corresponding to a paired difference of +0.116 (95% CI, +0.035 to +0.181) relative to observed clinician reward. However, the effective sample size was 51.5 events, indicating limited overlap between the learned and estimated historical policies. ConclusionsThis study demonstrates the feasibility of constructing and evaluating an offline contextual-bandit framework for lung donor-recipient matching using registry data. The estimated policy value was higher than the observed historical reward, but this finding is exploratory because candidate sets were reconstructed and off-policy evaluation had limited overlap. Verified offer sets, time-stamped candidate data, explicit compatibility constraints, and prospective evaluation are required before clinical use.

Estimating the correlation between subject-specific rates of change is a common objective in longitudinal studies. A natural approach is to fit ordinary least squares regressions per subject and correlate the resulting slope estimates directly. Under complete data, this slopes-first estimator is algebraically equivalent to a contrast-based plug-in estimator derived from the joint variance-covariance matrix of the two signals across time. Under independent dropout, however, the slopes-first approach implicitly applies heterogeneous weights across missing-data patterns, introducing bias whenever the cross-covariance between the two signals is not constant over time. We formalise this result and propose an unbiased alternative that estimates the joint covariance matrix before applying the contrast. We further extend our results to include covariates as well. Simulations confirm that the slopes-first estimator is biased under non-constant cross-covariance while the proposed alternative remains unbiased. Further, we show in a real-world example that the two estimators can lead to different scientific conclusions and recommend that researchers working with longitudinal data inspect the missing data pattern to make an informed decision about the estimator of choice.

Longitudinal non-Gaussian data are common in biomedical and microbiome studies, where responses may not be normally distributed. However, many existing dynamic graphical models are developed for Gaussian responses, while many non-Gaussian graphical models focus on cross-sectional dependence and do not explicitly estimate temporal transition effects. We propose a latent Gaussian dynamic graphical model for longitudinal non-Gaussian responses, including binary, count, and compositional outcomes. The proposed framework jointly estimates the contemporaneous conditional dependence network and the lag-one temporal dependence structure. We further develop a penalized Monte Carlo expectation maximization procedure for estimation. Simulation studies across multiple graph structures and outcome types show that the proposed method can recover both dependence structures, with more stable recovery for the lag-one temporal dependence structure and more structure-dependent recovery for the contemporaneous conditional dependence network. We apply the proposed method to longitudinal data from the Integrative Human Microbiome Project (HMP2) using a multinomial observation model. The results reveal persistent lag-one temporal dependence and contemporaneous conditional dependence among dominant microbial taxa, demonstrating the utility of the proposed framework for dynamic network analysis of longitudinal non-Gaussian biomedical data.

Biomedical Machine Learning (ML) is increasingly evaluated through Independent and Identically Distributed (IID) benchmarks, even when deployment involves irregular observation, missingness, distribution shift, and high failure costs. This position paper argues for a stability-oriented reporting standard for safety-relevant biomedical ML. The central requirement is not a specific model class or optimizer, but auditable evidence: performance should be reported under clinically plausible perturbations such as thinning, timestamp jitter, bursty missingness, channel dropout, and Signal-to-Noise Ratio (SNR) shift. When authors invoke training-stability diagnostics, those diagnostics should be hazard-linked: their assumptions, operational meaning, and relationship to stress-test degradation should be stated explicitly. We use controlled mechanism probes to illustrate two hazards that IID scores can hide: observation-process dependence under irregular sampling, and perturbation sensitivity associated with high-curvature training regimes such as Edge of Stability (EoS). We do not propose Neural ODEs, continuous-time models, Sharpness-Aware Minimization (SAM), or EoS metrics as default prescriptions. Instead, we propose a minimum reporting standard that makes robustness claims testable and supports escalation to more complex modelling or optimization only when stress tests justify the cost.

In systematic review, the measures reported by the primary studies may not always be the measures sought by the systematic review investigators. Motivated by an systematic review of interventions to reduce meat consumption, we discuss here methods to infer relative risks and confidence intervals from randomized controlled trials that reported differences, possibly in the presence of repeated measures or multiple interventions.

Background: Cost-effectiveness analyses (CEA) built on simulation models are usually delivered through project-specific scripts. Each new model, and often each new question, rebuilds the same scaffolding: parameter management, calibration, base-case analysis, sensitivity analysis and export. This is costly, and is where traceability between an estimate and the model version producing it is lost. We present THALASSA, a web application providing this infrastructure once, for any model satisfying a small programming contract. Methods: A simulation model is connected by one file implementing a fixed contract of R functions that determine its behaviour. THALASSA exposes the complete CEA workflow in a browser: a generated overview of the connected model, parameter editing, calibration, base-case comparison with the efficiency curve, and deterministic and probabilistic sensitivity analysis at a willingness-to-pay threshold, each with timestamped file export carrying the parameter values it was run with. Every session is bound to a repository version, so results carry their generating model version. Results: We demonstrate THALASSA on two public models of contrasting difficulty: a simple tutorial model and a generic cancer screening model with two competing natural-history pathways and a 21-dimensional calibration deliberately built to be more complex. Both ran on the same unmodified setup, implementing the same contract, with no code changes. Conclusions: Separating analysis infrastructure from model logic makes the CEA workflow reusable across projects and reproducible by construction. A modelling team can move from a versioned repository to audit-ready, exportable CEA outputs without writing analysis code. The integration effort does not grow with model size and one installation serves unrelated projects.

Understanding the causal pathways linking the gut microbiome to downstream biomarkers and clinical outcomes is central to microbiome research. When the microbiome acts as the exposure, the high dimensionality and compositional structure of the data introduce unique methodological challenges for mediation analysis. To address these challenges, we propose a latent variable mediation framework that captures variation in microbial composition through a microbial balance, defined as the log-ratio between two unknown subsets of taxa. This balance serves as a latent scalar exposure and simplifies the estimation of the overall indirect effect at the community level, while simultaneously identifying specific taxa that contribute to the overall direct and indirect effects. We describe the model's estimation and inference, and illustrate the method using data from the Multiethnic Cohort--Adiposity Phenotype Study, where we examine the mediation pathway from the gut microbiome through lipopolysaccharide-binding protein, a marker of metabolic endotoxemia, to percent liver fat.

Researchers must often rely on observational data to draw conclusions about the causal relationship between an exposure and an outcome. This is common in environmental epidemiology research, where the assignment to an exposure can be unethical or impossible. Consider a setting in which we seek to estimate excess mortality during heat waves, but are unable to randomize individuals to moderate or extreme heat. We propose to follow a causal framework for non-randomized time series that uses what we call historical time series matching to construct hypothetical randomized experiments, in which background covariates are balanced. We applied this framework and assessed the following causal question: Are daily mortality counts among elderly Americans higher during hot versus warm summer time periods? We addressed this causal question by drawing frequentist, Fisherian, and Bayesian inferences for five American cities. We estimated the average treatment effect of heat on three-day mortality with 95% Fisherian intervals to be 0.5 [-3.6, 4.4] deaths in Pittsburgh, 3.7 [0.7, 6.6] deaths in Seattle, 12.5 [6.1, 18.9] deaths in Chicago, 8.5 [1.4, 15.6] deaths in Los Angeles, and 20.3 [11.2, 29.5] deaths in New York.

Somatic mutations in sequenced tumors are often detected with imperfect sensitivity. They can be missed because tumor heterogeneity results in some mutations being absent in a sample at hand or present at very low frequencies which, combined with technical limitations like insufficient sequencing depth, reduces their detectability. It becomes apparent with multi-regional sequencing that some mutations confined to specific areas can be overlooked in single-region analyses. Adjusting for such sensitivity becomes critical when mutation frequencies between groups of subjects are compared, such as responders and non-responders to a given therapy. Otherwise, the estimate of the mutation prevalence, or of the difference in prevalences between the groups, can be biased, leading to inflated type I error. Here, we propose a weighted approach to adjust for imperfect sensitivity of the mutation detection and compared it to the unweighted approach. We assume that a surrogate of mutation sensitivity is observed only among patients with called mutation, and we derive its unconditional mean. Such adjustment leads to less biased estimates of mutational frequencies and better control of the type 1 error.

Accurate prediction of individual medical outcomes is essential for optimizing treatment allocation amid rising costs, coverage denials, and limited clinical resources. Traditional predictive models, including linear regression and neural networks, rely on average effects and cannot tailor predictions to the specific circumstances of individual cases. We present relevance-based prediction (RBP), a model-free and task-specific method that predicts outcomes as weighted averages of observed cases, with weights determined by a rigorously defined and theoretically justified measure of relevance. Unlike model-based methods that rely on fixed calibrated parameters, relevance-based prediction revisits the original data for each prediction and customizes both the cases and variables used. We illustrate RBP by applying it to predict opioid treatment outcomes and demonstrate that it provides case-specific insights unavailable from conventional models, including how each prior case informs a prediction, how each variable affects each predictions reliability and value, and how reliable each prediction is before it is made. These individualized insights may prevent misleading average-based decisions and reduce harmful or suboptimal treatment.

Mechanistic models of infectious disease transmission provide a framework for understanding and predicting disease dynamics, but statistical inference is often computationally prohibitive. This challenge is especially pronounced for models with high-dimensional structures and many unknown parameters, where trajectory matching may require thousands of evaluations of an expensive likelihood function. We developed a surrogate-assisted optimization framework that uses Bayesian optimization to reduce the number of evaluations of the underlying mechanistic model required for parameter inference. We evaluated the approach using simulated fitting problems based on age-structured susceptible-exposed-infected-recovered models with 4, 8, and 16 age classes and a 210-compartment age-structured rotavirus transmission model. Applying surrogate-assisted optimization substantially reduced the number of likelihood evaluations required to identify good parameter estimates, with larger improvements as model dimensionality increased. In the rotavirus model, the surrogate Bayesian optimization approach reached the benchmark after approximately 0.7 hours, compared with approximately 9.4 hours using traditional Nelder-Mead optimization. These results demonstrate that surrogate-assisted optimization can substantially reduce the computational burden of trajectory matching for complex mechanistic transmission models, particularly for higher-dimensional inference problems.

Doubly robust survival estimators require nuisance models that the available information can support. We examine design and estimator dependence in TCGA breast cancer data. Treatment indicators disagree for 337 of 1076 patients; 67.4% have unobserved five-year binary outcomes; and receptor imputation expands the triple-negative subgroup from 156 directly measured patients to 226. We compare chemotherapy initiation by day 90 with non-initiation, using adjusted restricted mean survival time (RMST) over the next 1095 days. Weighting-only estimates are 55 days (95% bootstrap interval: -24, 133) in triple- negative disease (16 deaths) and 38 days (12, 72) in the whole cohort (72 deaths). A Cox outcome-model benchmark yields 108 and 35 days, respectively, demonstrating material model dependence in the smaller population. A parameter-to-information ratio summarizes the burden on complete-outcome regressions within cross-fitting. Its simulation-calibrated warning level is 0.31 (0.28 - 0.35) for the original regime; additional simulations examine covariate-dependent censoring, delayed benefit and nonlinear misspecification. The ratio does not certify estimator validity or replace weight diagnostics. The contribution is a reproducible workflow that separates design, nuisance-model support and estimator stability before interpreting observational survival contrasts.

Motivation: Horizontal pleiotropy remains a major challenge for Mendelian randomization (MR), particularly as increasingly large genome-wide association studies identify many genetic instruments for complex traits. Existing pleiotropy-robust methods rely on assumptions regarding the proportion or distribution of valid instruments that may be violated when pleiotropy is widespread. Results: We developed CaPri, a Bayesian MR framework that jointly estimates causal and pleiotropic effects. CaPri allows a potentially arbitrary proportion of instruments to contribute to a shared directional pleiotropic component and estimates the posterior probability and magnitude of the causal effect, instrument-specific posterior inclusion probabilities for pleiotropy, and the shared pleiotropic effect. Across simulations, CaPri maintained appropriate false positive control and accurate causal effect estimation regardless of the proportion of pleiotropic instruments, while accurately identifying pleiotropic instruments and estimating their shared effect. Performance remained robust across varying causal and pleiotropic effect sizes and exposure-outcome confounding, although sufficiently large outcome GWAS sample sizes were important for stable inference and limited prior dependence. In real data analyses, CaPri reported low posterior probabilities of causality across negative controls and strong support for established positive control relationships. Application to alcohol consumption and coronary artery disease further demonstrated its ability to identify potentially meaningful pleiotropic instruments. CaPri provides a complementary framework for strengthening causal inference and triangulation across MR methods with different identifying assumptions. Availability and implementation: CaPri is available at https://github.com/OuRen0104/CaPri.

Background: Stepped-wedge cluster randomized trials (SW-CRTs) are increasingly used when interventions need to be implemented sequentially across clusters. Despite extensive developments in the design and analysis of superiority SW-CRTs, methods for planning noninferiority SW-CRTs remain limited. Methods: We propose swNIsim, a simulation-based framework for power calculation in noninferiority SW-CRTs using generalized linear mixed-effects models (GLMMs). The proposed approach accommodates both continuous and binary outcomes and supports practical design features, including unequal numbers of clusters across sequences, varying cluster-period sizes, delayed intervention effects, and analyses with or without adjustment for secular time effects. Noninferiority margins can be specified using several commonly applied strategies. Results: Through extensive simulation studies, we examined the impact of key design parameters on statistical power, including the number of periods, the number of clusters per sequence, the number of individuals per cluster-period, and the degree of between-cluster variability. For a fixed total sample size, designs with more time periods consistently achieved higher power. Increasing the number of clusters per sequence produced modest gains in power, whereas greater between-cluster variability reduced power. Type I error rates were well controlled across the scenarios examined. Balanced allocation of clusters across sequences was more efficient than unbalanced allocation schemes. Two real-world SW-CRT examples were used to illustrate application of the method. Conclusions: The proposed swNIsim provides a flexible and practical approach for power calculation in noninferiority SW-CRTs. By accommodating a wide range of design configurations and outcome types, and through implementation as a freely available R package, swNIsim may facilitate the planning and design of noninferiority stepped-wedge trials.

The human T cell repertoire is generated through recombination of a multitude of T cell receptor gene segments, yielding a complex array of T cell clones. T cell receptor beta (TRB) V gene segment defined, T cell clonal frequency when organized in correspondence with the respective V gene segment positions on the TRB genomic locus yields a periodic, undulating curve in the spatial domain. Using the genomic distance from the TRB-D1 segment to the TRB-V1-29 segments, spectral analysis was performed utilizing Lomb-Scargle periodogram (modified Fourier analysis) to obtain spectral power curves quantifying the TRB V clonal frequencies from 12 allogeneic stem cell transplant donors (baseline) and recipients (<100 days or b/w 100-365) using a variety of analytic software. Patients either underwent HLA matched related bone marrow or HLA matched unrelated blood stem cell transplantation, using PT-CY for the former and ATG for the latter. Spectral Power curves revealed dominant spectral peaks at wavelengths ranging from 113-335 millicycles/kb in the 12 donors, with consistent frequency domain spectral patterns, supporting similar use of V segments across healthy individuals. PT-CY recipients had power spectra closely corresponding to the respective donors, however the ATG recipients had relatively dispersed spectra, with spectral centroid shifted towards higher frequencies compared to donors and a reciprocal decline in the low frequency indices. Spectral power was concentrated in the <3 kb and 3-12 kb wavelengths in both groups, consistent with the periodicity observed in the relative V gene segment usage across the repertoire. The analyses reported here demonstrate that the healthy SCT donors have a spectral signature occupying short to intermediate wavelegnths in the frequency domain, PT CY recipients recover with donor-like spectra on the average, however ATG recipients tend to shift towards higher frequencies. These findings are consistent with a normal organized distribution of TRB V segment usage in healthy individuals, with PT-CY preserving this organization more consistently than ATG. Spectral (Fourier) analysis of TRB sequencing data provides a repertoire wide summary of T cell clonal distribution.

Background: Excess mortality estimates depend on an expected baseline, an unobservable quantity typically derived from models with prescribed functional forms and short reference periods that exclude exceptional events. Singular Spectrum Analysis (SSA) offers a nonparametric alternative but requires subjective component selection and is sensitive to outliers. We combined SSA with automated component selection and outlier interpolation to estimate long-term baselines and compile a catalog of excess mortality events (EMEs). Methods: We analyzed weekly all-cause mortality (2015-2024) for six metropolitan areas with contrasting climatic, demographic, and epidemiological contexts: Brussels, Antwerp, and Liege (Belgium), and Athens, Thessaloniki, and Larisa (Greece). SSA was applied to standardized log-mortality rates; components were selected using Monte Carlo SSA, and COVID-19 periods were interpolated iteratively. EMEs were defined as weeks exceeding a quasi-Poisson two-z-score threshold, then screened for coincidence with heat waves, cold waves, COVID-19 periods, and influenza epidemics. Results: Baselines were consistent within countries: annual cycles in Belgium, with an additional significant semiannual cycle in Greece reflecting summer mortality. We cataloged 170 EMEs spanning 456 weeks and 14,210 excess deaths. Pandemic weeks accounted for 77% of excess mortality, and 96 weeks coincided with temperature extremes; 56 EMEs, representing 3.7% of excess mortality, remained unmatched. Conclusions: SSA yields reproducible, low-assumption baselines that capture country-specific seasonality. Because flexible baselines can absorb recurrent seasonal mortality, component selection is a substantive choice that should be reported explicitly, together with its sensitivity.

In recent years, the issue of reproducibility in scientific studies has been at the forefront of scientific discussion. Many of these issues arise due to poor employment of statistical methods, making it crucial to address these knowledge gaps for frequently used methods. Elastic net regression (ENET) is one such algorithm used frequently in biology and public health due to its ability to select important predictors (also known as "features") from large data sets. However, ENET struggles with consistent feature selection and hyperparameter choice in these applications due to low observation counts and noisy data. This limitation hinders ENET's effectiveness in biological research due to low model reproducibility and predictor interpretability. We propose the Instability Quotient Criterion (IQC) as a metric to quantify the consistency of selected features for feature selection models such as ENET, with a specific emphasis on biological applications. IQC leverages ensemble modeling and Shannon's entropy in order to determine the stability of selected features derived from ENET models. We demonstrate IQC's robustness via simulation studies, showing that IQC is capable of detecting significant decreases in model stability with lower observation counts, higher dimensionality, and higher proportions of noisy signal. We additionally utilize simulation to provide a naive threshold for low versus high stability. Finally, we applied IQC to a case study to illustrate how our method can be integrated into a larger analysis pipeline when looking at microarray gene expression data. IQC indicated stability issues for differentially-expressed genes when classifying acute leukemia subtypes in spite of high classification success. Lastly, we discuss how model instability impedes interpretation of biologically significant genes for discerning cancer subtypes. Overall, we present IQC as a robust evaluation metric that can easily be incorporated into high-dimensional biological analyses.

Next-generation sequencing data is compositional in nature and ignoring this fact can significantly impact downstream data analyses. This is especially true in co-occurrence (correlation) networks because compositional data is prone to many false negative correlations. The centered log-ration (clr) transformation is a commonly used compositional data analysis (CoDA) approach that converts compositionally constrained data into real number space. The clr transformed values can subsequently be used in machine learning and statistical analyses such as Spearman or Pearson correlations used to generate co-occurrence networks. Here we show that the necessary zero substitution performed prior to clr transformation can create numerous multicollinearity artefacts. We also show that a newer CoDA approach called the robust centered log-ratio (rclr) eliminates these artefacts in correlation-based network analysis.